Long‐term community shifts driven by local extinction of an iconic foundation species following an extreme marine heatwave

Abstract Gradual ocean warming combined with stronger marine heatwaves (MHWs) can reduce abundances of foundation species that control community structures, biodiversity, and ecosystem functioning. However, few studies have documented long‐term succession trajectories following the more extreme events that cause localized extinctions of foundation species. Here, we documented long‐term successional changes to marine benthic communities in Pile Bay, New Zealand, following the Tasman 2017/18 MHW, which caused localized extinctions of dominant southern bull kelp (Durvillaea sp.). Six years on, multiscale annual and seasonal surveys show no sign of Durvillaea recolonization. Instead, the invasive annual kelp (Undaria pinnatifida), rapidly colonized areas previously dominated by Durvillaea, followed by large changes to the understory community, as Durvillaea holdfasts and encrusting coralline algae were replaced by coralline turf. Between 3 and 6 years after the total loss of Durvillaea, smaller native fucoids colonized in high densities. Although Undaria initially colonized plots throughout Durvillaea's tidal range, later in the succession Undaria only retained dominance in the lower intertidal zone and only in spring. Ultimately, the tidal zone was slowly replaced by alternative foundation species, composed of different canopy‐forming brown seaweeds that dominated different intertidal elevations, resulting in a net increase in canopy and understory diversity. This study is a rare example of long‐term effects following an extreme MHW that caused extinctions of a locally dominant canopy‐former, but these events and their associated dramatic changes to community structures and biodiversity are expected to become increasingly common as MHWs continue to increase in strength, frequency, and duration.


| INTRODUC TI ON
Recent studies have documented intensification of marine heatwaves (MHWs; Frölicher et al., 2018;Oliver et al., 2019Oliver et al., , 2021Thoral et al., 2022) and their impacts on marine organisms, including competitively dominant foundation species that control community structures (Doney et al., 2012;Smale et al., 2019). Impacts from MHWs on the ecological performance of foundation species has typically been reported on canopy-forming seagrass (Arias-Ortiz et al., 2018;Strydom et al., 2020) and different large brown seaweed (Arafeh-Dalmau et al., 2019;Filbee-Dexter et al., 2020;Rogers-Bennett & Catton, 2019;Straub et al., 2019;Tait et al., 2021;Wernberg et al., 2016). During these MHWs, temperatures may exceed physiological thermal optimums resulting in decreased ecological performances and lowered species abundances (Filbee-Dexter et al., 2020;Rogers-Bennett & Catton, 2019;Smith et al., 2023;Spiecker & Menge, 2022;Suryan et al., 2021;Tait et al., 2021;Whalen et al., 2023). However, a few rare case studies have documented much more conspicuous ecological effects, like when extreme MHWs cause regional extinction of the competitively dominant foundation species Wernberg et al., 2016). It is likely that localized extinctions of competitively dominant foundation species will become more common because MHWs are predicted to become stronger, longer, and more frequent . For example, over the last 40 years, most temperate coastal ecosystems, which often are dominated by large brown seaweeds, have experienced strong increases in MHWs (Thoral et al., 2022). More specifically, similar trends have been observed in coastal New Zealand (Montie et al., 2023), a bioregion characterized by a long temperate coastline, high biodiversity, and unique marine communities with many endemic species (Costello et al., 2010). It is therefore of principal importance to document how coastal ecosystems may change in the future following the total loss of competitively dominant foundation species.
The wave-exposed intertidal and shallow subtidal coastline on New Zealand's South Island is often dominated by the native foundation species, Durvillaea antarctica, D. willana, and D. poha (the two latter species are also endemic), more generally referred to as southern bull kelp. Southern bull kelp can reach 10 m in length, live for 7-8 years, build canopy-dominated seascapes, control community structure and productivity, and are therefore considered iconic species that epitomize New Zealand's rugged coastline and kelp-associated fisheries of abalone, butterfish, and crayfish (Fraser et al., 2009;Hay, 1994;Taylor & Schiel, 2005;Velásquez et al., 2020). During the summer of 2017/18, the east coast of the South Island experienced extreme air and sea surface temperatures, corresponding to the strongest MHW on record (Salinger et al., 2019;. This extreme climatic event began on 14 November 2017 and lasted 147 days until 9 April 2018 (Salinger et al., 2019). It resulted in rapid extinction of bull kelp at our study area, Pile Bay, Lyttelton Harbour, and at other sites within and around Lyttelton Harbour . This extinction event was followed by immediate colonization of the invasive kelp Undaria pinnatifida (hereafter Undaria), one of the world's most widely distributed marine invaders that has inhabited Lyttelton Harbour since the early 1990s (Epstein & Smale, 2017;South et al., 2017;Thomsen et al., 2018. Here, we provide a rare case study of the long-term (6 years) successional trajectories after complete loss of southern bull kelpthe locally competitively dominant foundation species-and colonization by different native and invasive foundation species-across spatiotemporal scales. The few existing case studies suggest that replacement of lost foundation species with alternative foundation species is a slow process Wernberg et al., 2016). Considering bull kelp were extirpated on all reefs in Lyttelton Harbour and nearby , recolonization would rely on drift of separate male and female fronds during a short reproductive period from autumn to early spring (Fraser et al., 2018;Taylor & Schiel, 2003;Velásquez et al., 2020). We therefore hypothesized that bull kelp would not recolonize Pile Bay and that species that typically co-occur with bull kelp, like encrusting algae, would also decrease in abundance . We also hypothesized that Undaria would continue to benefit from reduced competition from bull kelp and that space previously dominated by bull kelp eventually would be colonized by other seaweed that are common on this particular reef, like coralline turf and other native brown seaweeds .

| Regional context
Surveys were conducted along the rocky intertidal zone in Pile Bay (−43.615, 172.765), Lyttelton Harbour, Banks Peninsula. Pile Bay is a less wave-exposed reef compared to typical bull kelp-dominated reefs. The reef is composed of volcanic rock (Sewell, 1985) and seasonal SST ranges between 9 and 20°C. However, during the Tasman 2017/18 MHW, SST increased to more than 23°C (Salinger et al., 2019;. This was the strongest MHW event recorded over the 38 years of existing satellite data, with maximum and cumulative intensities two times greater than any recorded event since 1980 . The focal foundation species of this study, southern bull kelp, are adapted to wave-exposed cool, clear, nutrient-rich waters (Hay, 1994;Velásquez et al., 2020), potentially making them vulnerable to anthropogenic stressors like warming or decreased water clarity. The pre-MHW bull kelp population in Pile Bay all had wide blades, a yellow-brownish tint, and were relatively short (<5 m), suggesting that they were dominated by the endemic D. poha (and not D. antarctica;Fraser et al., 2012;Velásquez et al., 2020).

| Changes to dominant seaweed on the reef scale (drone images)
To test whether the Tasman MHW affected the general distribution   patterns of intertidal canopy-forming seaweed in Pile Bay, large-scale   drone surveys were conducted before the MHW in August 2017 with an Advanced Phantom 3 equipped with a 12 MP, 1/2.3″ CMOS sensor and FOV 94° 20 mm lens with an image size of 4000 × 3000 and ground sample distance of 0.43 cm/pixel. Surveys were repeated annually until September 2022 with a Mavic Mini 2 equipped with a 12 MP, 1/2.3″ CMOS sensor and FOV 84° 24 mm lens with an image size of 4000 × 3000 and ground sample distance of 0.36 cm/pixel. Geotagged RGB images were taken at an altitude of 10 m during spring tides (0.1-0.3 m above lowest astronomical tide [LAT]), where each image covered c. 95 m 2 . For each image, percent cover of large canopy-forming algal species was visually estimated using a superimposed grid of 100 cells and separated into the low, mid, and high intertidal zones, based on cover of water and bare substrate (cover of water between 10% and 40% = low zone, <10% water and <50% bare substrate = mid zone and 0% water and >50% bare substrate = high zone). Changes to community structure and percent cover of the dominant canopy-forming seaweed were analyzed with two-way permutational multivariate analyses of variance (PERMANOVA) using Bray-Curtis dissimilarity coefficients for multivariate community data (see Appendix S1 for results) and Euclidian distance for univariate taxa responses. The two fixed factors included year ( pre-MHW, 2018( 1-year, 2019( 2-year, 2020( 3-year, 2021( 4-year, and 2022 5-year post-MHW) and elevation (low, medium, and high). Data could not be transformed to achieve homogeneous variances, so results should be interpreted with caution. We therefore adjusted alpha to .01 and acknowledge that significant effects may alternatively have been caused by heterogeneous variances (Underwood et al., 1997).
However, because we were interested in evaluating interaction effects and assessing relative importance of crossed test factors (Levine & Hullett, 2002), biased factorial PERMANOVAs were preferred over nonparametric ranked tests that do not address interaction terms or relative importance of test factors (but still require similar distributions). Furthermore, PERMANOVA approaches are generally robust to minor heteroscedasticity for balanced designs (Quinn & Keough, 2002;Underwood et al., 1997). A similarity of percentages breakdown (SIMPER) was used to assess the contribution of individual species to community dissimilarity between years and elevations.

| Common species in pre-heatwave bull kelp forest (understory photographs)
Before the MHW in August 2017, six 1 m 2 quadrats with 100% cover of bull kelp canopy were photographed at two bull kelp-dominated sections of the reef (see Figure 1, i.e., a total of 12 quadrates). Before a photograph was taken, the large bull kelp fronds were pushed aside so that abundances of understory species could also be quantified.
Within each photograph, percent cover of sessile organisms and counts of mobile organisms were recorded to their lowest taxonomic level (genus or species). Photographs were taken from random and scattered bull kelp plots covering the low-to-mid elevation zone on both reef sections, but plot-specific elevation levels could not be determined from image analyses (and results are therefore only reported for the combined low-mid elevation level). Furthermore, because pre-MHW data were estimated from these photographs, some cryptic fauna like very small chitons, limpets, and snails, which can inhabit narrow cracks and crevices, may be underestimated compared with follow-up post-MHW in situ surveys (see below).

| Impact of Undaria, elevation and season on succession (permanent plots)
To test how tidal elevation and invasive Undaria affected succession trajectories following the bull kelp extinction event, 0.25 m 2 (0.5 × 0.5 m) permanent plots were established in habitats previously dominated by bull kelp on the two reef sections photographed in 2017.
In March 2018, 12 permanent plots were marked on reef section 1 and 16 plots on section 2, with another 4 plots (that were not established initially due to bad weather, poor tides, and logistic problems) marked on reef section 2 in September 2020, to ensure the two reef sections had similar number of plots and sampling effort ( Figure 1a).
The elevation level of each plot above LAT was measured to nearest cm with an Abney level and reclassified to low (<0.7 m above LAT) or mid (>0.7 m above LAT) intertidal zones. In half of the permanent plots (eight per site), an Undaria removal treatment was applied where colonizing Undaria were regularly removed before sampling, and during sampling if recruits were found, to test whether this invasive spe- Changes to community structure were visualized with nMDS, by combining data from the pre-MHW photographs and the post-MHW permanent plot surveys. Three-factorial PERMANOVAs were done on the post-MHW data only (because they were sampled differently than the pre-MHW data) using annual spring data on multivariate community structure (using Bray-Curtis dissimilarity coefficients) and univariate responses (using Euclidian distances). Univariate responses included the abundances and taxonomic richness of all sessile (percent cover) and all mobile (counts) species, and key habitatforming species (percent cover), that is, Undaria, Hormosira banksii, Cystophora torulosa, C. scalaris, coralline turf, and encrusting coralline algae. Year (1-, 2-, 3-, 4-, and 5-year post-MHW) was considered a fixed factor and the factors Undaria removal (control and removal) and elevation (low and mid) were nested within a random plot factor. We supplemented the annual tests with repeated measures four-way PERMANOVAs on the seasonal data (autumn and spring 2020-21 and 2022-23 using the same parameters as the three-way PERMANOVA). However, the seasonal tests were less important for our hypothesis related to long-term changes and were therefore only reported in online supplements. Again, data could not be transformed to variance homogeneity, but parametric PERMANOVA was (again) preferred over nonparametric ranked approaches for reasons detailed in Section 2.1. Significant test results were followed up with post hoc t-test comparisons to identify significant treatments.

| Habitat interactions (permanent plots)
To test whether habitat interactions in the post-MHW community were affected by elevation and season "attachment surveys" were done in the 32 permanent plots, using a quadrat of 0.25 × 0.25 m (0.0625 m 2 ) over 2 years from 2021 to 2023 (two times in spring and two times in autumn). For this survey, the small quadrat was placed in a random subsection of the larger 0.25 m 2 permanent plot.
Cnidarians and bryozoans remained low in abundance ( Figure S2).
Finally, mobile mollusks substantially increased after the MHW reaching highest levels in spring 2022-23 (28 ± 2 SE at low elevation and 20 ± 4 SE at mid elevation, Figure S3).

F I G U R E 2
Average abundance (percent cover ± standard error) of four key habitat-forming fucoids (Cystophora sp., including C. torulosa and C. scalaris, Durvillaea sp. and H. banksii) and kelp (Undaria) detected from reef scale drone imagery. Drone surveys were conducted annually between 2017 and 2022 during spring, when Undaria cover was highest. Average percent cover was calculated based on transects along the low (blue), mid (green), and high (red) intertidal zones. Data were not captured in the high intertidal zone in 2017 or 2019. The gray shading represents pre-MHW coverages.
F I G U R E 3 NMDS plot showing changes in community structure at sites previously dominated by bull kelp using pre-MHW (2017) data from photographs of 1 m 2 quadrats (n = 12) in the intertidal (black squares) and post-MHW data from permanent 0.25 m 2 plots (n = 32) separated into low (circles) and mid (triangles) intertidal elevations. The pre-MHW data have a gray-shaded background to highlight their extraction from photographs (where specific elevation levels could not be determined), whereas post-MHW data were estimated from visual surveys of permanent plots and highprecision elevation measurements.

| Taxonomic richness and abundances of mobile and sessile species
The three main test factors explained in concert 69% (sessile richness), 56% (mobile richness), 72% (sessile abundance), and 68% (mobile abundance) of the total data variability. Overall, there was relatively high interplot variability (i.e., Plot (Undaria removal × elevation) tests were generally significant). We found no effect of elevation or Undaria removal (p > .01) and only a single significant interaction that explained little of data variability (Sessile abundance: p year × elevation = .002, 5% of SS, Table 2D). By contrast, year was, as expected, significant across responses and explained a high proportion of the data variability (p sessile richness = .001, 51% SS, p mobile richness = .001, 17% SS, p mobile abundance = .001, 35%, p sessile abundance = .001, 28%, Table 2A

| Abundances of key species
After the MHW, percent cover in the permanent plots was dom-  (Table 3). Again, we found high interplot variability across tests (all Plot (Undaria removal × elevation) tests were significant) but only F I G U R E 4 Average abundance (percent cover ± standard error) of key habitat-forming fucoids (Durvillaea sp., C. torulosa, C. scalaris, and H. banksii) and kelp (Undaria), understory coralline turf (C. officinalis) and encrusting coralline algae, in photo quadrats before (2017, n = 12 1 m 2 random plots, black points, gray shading, no elevation measurements) and after (n = 32 permanently marked 0.25 m 2 plots, 2018-2022/23) the MHW, where post-MHW data were separated into low (blue) and mid (green) intertidal elevations during spring (circles) and autumn (triangles). The abundance of Undaria only includes the data from the control plots where Undaria was not experimentally removed.

| Habitat interactions (permanent plots)
We found more attachment interactions in the low compared to mid zone and in autumn compared to spring (Figure 6a), whereas the analyses on the number of linkages revealed more complicated results ( Figure 6b). The three test factors explained in concert 7% and 5% of the data variability, for number of attachment and linkages, respectively (Table 4). We found significant interactions for year × elevation for attachments (p = .011, 1% of SS) and year × season for linkages (p = .001, 2%). Furthermore, year was only significant for linkages (p = .006, 0.6%) whereas linkages and attachments were both affected by season (p linkages = .001, 0.7%, p attachments = .001, 5%) and elevation (p attachments = .001, 0.6%, p linkages = .001, 1%) (Table 4).
Most attachment chains involved only two individuals, but a few had three or more species, for example, the three-chain interaction with H. banksii-attached-to-C. maschalocarpum-attached-to-C. officinalis-attached-to-encrusting coralline alga.

| DISCUSS ION
The rise in extreme MHWs across the globe, and their ecological impacts, is of growing concern, especially when these events cause the loss of habitat-forming foundation species (Arafeh-Dalmau et al., 2019;Straub et al., 2019;Tait et al., 2021;Wernberg, 2021;Wernberg et al., 2016). Here, we investigated community shifts after an extreme MHW resulted in the local extinction of an important foundation species. The community was surveyed over 6 years and impacts associated with annual, seasonal, and elevational test factors were quantified. We found that the extreme MHW caused long-lasting, possibly permanent, community changes to a coastal ecosystem after the dominant foundation species was lost, triggering a shift in understory species, and followed by high recruitment and replacement by alternative foundation species, that, in concert, likely will inhibit future recolonization of lost  (2018, 2019, 2020, 2021, and 2022) was fixed and factors Undaria removal (control and removal), and elevation (low and mid) were nested within the random factor plot (a repeated measure). Significant p-values are shown on bold (alpha < .01; variances were heterogenous).
Globally, research suggests that localized extinctions can result in slow recovery of the once-dominant canopy-former . However, here we found support for our hypothesis that bull kelp did not return to Pile Bay reefs during a sampling period of 6 years, and bull kelp recruits are yet to be found in this region (pers. obs. May 2023, the nearest surviving population remain >10 km away). While bull kelp has effective settlement of fertilized eggs over short distances (meters) and short time scales (hours), hardly any zygotes disperse longer distances (Fraser et al., 2018;Taylor & Schiel, 2003;Velásquez et al., 2020). Therefore, reef-toreef recolonization is more likely to arise from detached fertile male and female thalli arriving as intertidal beach cast at the exact same time and place (Velásquez et al., 2020). Furthermore, colonization by coralline turf, Undaria and smaller canopy-formers, can inhibit bull kelp recruitment because zygotes and juveniles settle and survive best on bare rock and encrusting alga (Schiel, 2019;Taylor & Schiel, 2005). In short, it is not surprising that bull kelp have not recolonized Pile Bay, even though this location is positioned well within D. poha latitudinal range (its northern range limit is 200 km north; Fraser et al., 2018Fraser et al., , 2012Taylor & Schiel, 2003;Velásquez et al., 2020). Similar severe and long-lasting impacts from a MHW have been documented in Western Australia, where lost foundation species have not recolonized 8 years after an extreme event caused c. 100 km range contractions of the kelp Ecklonia radiata and the fucoid Scytothalia dorycarpa Straub et al., 2019;Wernberg, 2021;Wernberg et al., 2016).
In our study, bull kelp was replaced by different seaweed and animals at slightly different elevations, which resulted in higher net biodiversity and abundance of mobile fauna and sessile understory seaweeds (minus those associated with bull kelp like encrusting coralline algae). Some of these species are "alternative foundation species," which, like bull kelp, provide habitat for plants and animals and alter abiotic conditions, for example, by reducing light and increasing relative humidity in the subcanopy (Elsberry & Bracken, 2021;Lilley & Schiel, 2006;. However, in contrast F I G U R E 6 (a) Average number of habitat attachments (± standard error, proxy for abundance) and (b) Average number of linkages (± standard error, proxy for richness of species involved in habitat interactions) within permanent plots during seasonal surveys in spring (red) and autumn (black) between 2021-22 and 2022-23.
to heavy and large-bladed bull kelp, canopies of the alternative foundation species do not "whiplash" the substrate preventing some individuals from inhabiting the substrate underneath (Stevens et al., 2002;Taylor & Schiel, 2010), partly because they are much smaller (0.5-2 m). Indeed, bull kelps are the largest fucoid seaweed worldwide, with fronds up to 10 m in length and biomass reaching 80 kg m −2 (Santelices et al., 1980;Taylor & Schiel, 2005), and the fronds create a whiplash effect under wavy conditions that maintain a unique understory of encrusting alga, certain epiphytes and mollusk, like abalone and chitons (Santelices et al., 1980;Schiel, 2019;Stevens et al., 2002;Taylor & Schiel, 2010;. Furthermore, although the alternative foundation species benefitted after bull kelp was lost, endemic D. poha has a narrow latitudinal range (Fraser et al., 2012) and appears to be vulnerable to extreme MHWs , which will have cascading effects on diversity and abundance of the many culturally and economically important species they support, like fish and abalone, as well as the multitude of grazing invertebrates that reside in their large holdfasts (Fraser et al., 2011;Hay, 1977;Santelices, 1990;Schiel et al., 2018;Taylor & Schiel, 2005).
Intertidal habitats vary dramatically in environmental conditions across short vertical gradients, as desiccation and temperature fluctuations increase with shore height, but competition for space and predation by marine consumers decreases (Davenport & Davenport, 2005;Raffaelli & Hawkins, 1996;Schreider et al., 2003;Thomsen et al., 2016). However, our analysis showed relatively minor effects of elevation on new colonization by alternative foundation species, although a more diverse benthic understory community developed 3-4 years after the extinction event, with vertical (niche) partitioning as different subordinate species colonised slightly different elevation levels (that previously was dominated by bull kelp).
For example, coralline turf and Undaria immediately colonized and dominated in the mid-intertidal zone, whereas interspersed H. banksii and Cystophora sp. became the dominant habitat-formers in the mid-to-high intertidal zone, as found in past clearance experiments (Lilley & Schiel, 2006;Schiel & Lilley, 2011). We also found, as hypothesized, that Undaria continued to benefit from open spaces across elevations until 2019, whereafter bare space decreased (as more species colonized), and Undaria became less common in higher tidal elevations-eventually matching its typical vertical distribution near the low water neap tide mark (Forrest & Taylor, 2002;South et al., 2015). Indeed, Undaria was the only alternative foundation species that was significantly affected by elevation. To our knowledge, mechanisms driving Undaria's upper vertical limit have not yet been experimentally investigated but are likely caused by desiccation stress and perhaps also competition with other intertidal algae and invertebrate grazing (South et al., 2017).
We also showed that season significantly affected sessile taxonomic richness and mobile species abundance, as well as the abundance of four key habitat-forming species and understory corallines (see Appendix S1), probably because of cyclic changes in external environmental conditions like dissolved nutrient concentration, temperature, light availability, and wave action (Dayton, 1985;Lilley & Schiel, 2006;Schiel & Lilley, 2011). Indeed, large perennial brown alga, like kelp and fucoids, have complex life histories evolved to respond to these predicable changes (Clayton, 1988;Schiel & Foster, 2006).
For example, invasive Undaria is a winter annual species with a sporophyte that grows rapidly during winter and spring forming relatively short-lived but dense monospecific stands (South et al., 2017 impact on native species successional trajectories, probably because of the species pronounced seasonality, that allowed native species to colonize plots in spring and summer (when Undaria is sparse or absent). Finally, we showed that species-attachment interactions were more common in the low zone during autumn when plots have lower Undaria cover but higher cover of mixed H. banksii and Cystophora sp. canopies, and understory coralline turf, where the latter two taxa provide habitat for many epiphytes (like the brooding anemone Cricophorus nutrix) and mobile snails and whelks, like in . Unfortunately, we did not have habitatinteraction data from the pre-MHW bull kelp bed, but other research suggests that bull kelp typically facilitate encrusting alga (by whiplashing more competitively dominant seaweed), which subsequently provide habitat for a few red seaweed like Ballia spp., Gelidium microphyllum and various limpets and chitons through facilitation cascades (Schiel, 2019;Taylor & Schiel, 2005;. These specific bull kelp community interactions have largely been lost after the MHW and are instead replaced with different interaction networks now driven by Cystophora sp. and coralline turf . Furthermore, these types of attachment interactions highlight the importance of facilitation cascades in modifying and enhancing marine biodiversity (Altieri et al., 2007;Gribben et al., 2019;Thomsen et al., 2010).
In conclusion, we documented long-term successional community changes following the extirpation of a temperate marine foundation species, driven by an extreme MHW. Specifically, lost endemic and native bull kelp did not recover over 6 years but was instead replaced by invasive and native alternative foundation species that altered interaction networks and increased diversity of understory flora and fauna, like coralline turf, that can inhibit future recolonization of bull kelp. As MHWs are predicted to become longer and stronger, localized extinction events and similar conspicuous changes to communities and associated ecosystem functions and services (Smith et al., 2021) are expected to become increasingly common and will be a major threat to endemic cold-water species with narrow latitudinal ranges, like D. poha.

CO N FLI C T O F I NTER E S T S TATEM ENT
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are openly available in "figshare" at https://doi.org/10.6084/m9.figsh are.23519 679.v1.